Split-step Fourier methods for the Gross-Pitaevskii equation 
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We perform a systematic study of the accuracy of split-step Fourier transform methods for the 
time dependent Gross-Pitaevskii equation using symbolic calculation. Provided the most recent 
approximation for the wave function is always used in the nonlinear atom-atom interaction potential 
energy, every split-step algorithm we have tried has the same-order time stepping error for the Gross- 
Pitaevskii equation and the Schrodinger equation. 
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It is now firmly estabfished that the time-dependent 
and independent versions and muhi-component gener- 
alizations of the Gross-Pitaevskii equation (GPE) give 
a useful basic picture of most of the phenomena taking 
place in a dilute atomic Bose-Einstein condensate P, U 
At this writing there are hundreds of journal publica- 
tions involving direct numerical solutions of the GPE. 
The same equation, known as the nonlinear Schrodinger 
equation, is a central topic in nonlinear optics 3]. 

Curiously, many research groups have their own unique 
numerical methods, the methods are often ad-hoc, 
and the convergence properties are not always stated 
(known?) explicitly. Even if we restrict the discussion 
to the usual time dependent GPE, there is enough va- 
riety to thwart our attempts at a classification; but the 
Crank-Nicholson scheme and various split-step methods 
are common themes J3||. Our focus is on the analogs of 
the classic split-step Fast Fourier Transform method for 
the Schrodinger equation |^ . The idea is to split the GPE 
equation into two parts, one of which only refers to the 
momentum operator and the other only to the position 
operator. The wave function is then evolved alternatingly 
in momentum space and in real space. 

Nonetheless, the original split-operator method is 
formally based on the algebra of linear operators, and 
does not directly go over to the nonlinear GPE. Here we 
study the accuracy of split-operator methods for the GPE 
using power series expansions in the time step. While 
straightforward in principle, with increasing order such 
expansions rapidly becomes unwieldy, and are only man- 
ageable using symbolic calculation. Our surprise discov- 
ery is that there is a simple and computationally inex- 
pensive way to include the nonlinearity so that the split- 
operator method works for the GPE, and all the favorable 
properties of the original algorithm that have made it 
a huge success are preserved. 

We first recap the case of the time-dependent 
Schrodinger equation, for the time being in one spatial 
dimension. For brevity of the notation we use units such 



that the mass of the particle and the constant h both 
equal unity. The problem is to integrate 
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T = -^^, V = V{x). (2) 

A time step from t to t + h is carried out formally as 

^{x,t + h)^e-'''^^+^^^{x,t). (3) 

In position representation e^*''^ means multiplication by 
the function g-'^f^'^i^) ^ Similarly, in momentum represen- 
tation, after the Fourier transform J^[ip(x,t)] —^ ip{p,t), 
the kinetic-energy exponential multiplies the wave func- 
tion ip{p,t) by e~^'^P 1"^ . However, the operators T and 
V do not commute, so that the inequality e~*'''-^+^) ^ 
g-i-hv g-thT }^Qi([g g.j^^ ^jjg exponential of the sum of the 
two operators may be difficult to calculate. 

Split-step methods attempt to get past the obstacle of 
noncommuting operators by approximating 



(4) 



In the present discussion it does not matter what the 
linear operators are, so we refer to generic A and B. 
In practice the split is not useful unless there is an 
easy way to calculate each operator exponential on the 
right-hand side, but this too is immaterial in our formal 
development. Finally, we need not be dealing with time 
stepping, thus we write the scalar parameter as A. In 
general, though, we regard A as small in absolute value. 
The idea of split-operator methods is to pick the coeffi- 
cients ai and Pi so that the right-hand side of the split Q 
approximates the left-hand side to as high an order in A 
as possible. 

We seek split-operator methods using symbolic manip- 
ulation on Mathematica 6] . We first expand the exponen- 
tials, taking care not to inadvertently swap the noncom- 
muting operators A and B. Next, in order to compare the 
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TABLE I: Nonzero coefficients for minimal split-step 
methods with real coefficients for the orders of error 
0(X^),0{X'^), and 0{\^). The temporary notations are V = 
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left and right sides of Eq. Q), it is expedient to put the 
operators into a standard order, e.g., all operators B to 
the left and all operators A to the right. This introduces 
commutators of the operators, as in AB = BA + [A, B]. 
The difference between the left and right sides of Eq. Q 
is then arranged into a power series of the parameter A, 
where the coefficients contain commutators of the opera- 
tors and the constants a^. Pi. The requirement that the 
difference cancels order by order in A gives multivariate 
polynomial equations in ai and which Mathematica 
appears to solve easily and completely. 

The simplest nontrivial split-operator method found 
in this way is the original split-operator algorithm 5] 
with three exponentials. The general operator- algebra 
argument does not distinguish between the operators A 
and S, and the choices for the nonzero coefficients ai — 
ai = 1/2, /3i = 1 and /3i = = 1/2, a-i = \ wiU both 
do. We have, for instance, a split-operator representation 
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Time stepping as in Eq. (j3Jl is readily implemented using 
the Fast Fourier Transformation. The ensuing algorithm 
automatically preserves the norm of the wave function, 
which is an issue when standard differential equation 
solvers are applied to the time-dependent Schrodinger 
equation. The implementation of the exponential of the 
kinetic energy is a high-order spectral method, and the 
exponential with the potential energy is done exactly in 
principle. 

A five-exponential split exists with an error 0(A*), but 
the coefficients cti and /?i are not real and the algorithm 
is not absolutely stable. There is a whole family of six- 
exponential splits with a continuous- valued free parame- 
ter that all have an error 0{\^\ and the coefficients are 
real for a range of the values of the free parameter. It 



takes a seven-exponential split to gain another reduction 
of the error to ©(A^), and higher-order methods are also 
found without difficulty. However, for the present pur- 
pose we stop at ©(A^), and list in Table^the coefficients 
Q!i, Pi for all minimal (as few exponential factors as pos- 
sible for a given order) split-operator methods with real 
coefficients up to the order O(A^). The first publication 
known to us that gives these coefficients is Ref. p]. 

On the other hand, the Gross-Pitaevskii equation 
(GPE), or the nonlinear Schrodinger equation, reads in 
suitable units 
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where g (usually > 0) encapsulates the strength of the 
(usually repulsive) atom-atom interactions. The GPE 
is nonlinear and as such it does not fit the operator al- 
gebra framework we have used to derive split-operator 
methods, but the square of the wave function formally 
behaves like a potential energy. One is tempted to add 
it to the potential-energy term in the time stepper, as in 
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Nonetheless, as the wave function evolves in time, it 
is not clear what to insert for the wave function in the 
exponential. In fact, doing the step as in Eq. (0 drops an 
order in accuracy, and the error turns out to be 0{h^). 

In an attempt to find a better method we split the 
split-operator step explicitly, 
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At the stage when the square of the wave function is 
needed, there are already two versions of the wave func- 
tion available. The hope is to pick a linear combination 
of the two with the so far unknown coefficients cq and ci 
to regain the error 0{h^). 

In order to produce the "exact" result for comparisons 
we write 

h d /i2 52 

V'(x, t + h)^ i^{x, t) + - —^/j{x, t) + — Q^^{x, t) + .... 

(9) 

The nth time derivative of the wave function is obtained 
inductively, by taking the (n— l)th time derivative of the 
GPE and using the already existing expressions for the 
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derivatives up to the order n — 1 to ehminate all time 
derivatives from the right-hand side. This procedure is 
carried out using Mathematical treating real and imagi- 
nary parts of the wave function separately. A mechani- 
cal implementation is well advised, as the terms multiply 
rapidly with the order n; after complete expansion of 
the derivatives, the real part of the fourth time deriva- 
tive already has 168 linearly independent terms. The 
split-operator algorithm to be tested is implemented by 
expanding the exponentials into power series of the oper- 
ators, and acting the series-form operators on the initial 
wave function in the prescribed sequence. Such expan- 
sions also become extremely tedious when the order in 
h increases, and they are done using Mathematica. The 
result is that the method (jHJ gives an error 0{h^) if and 
only if the coefficients cq and ci satisfy cq = 0, |ci| = 1; 
in other words, if the most recent available wave function 
is used in glV'P- 

We have done a similar analysis for all split-operator 
methods with the coefficients listed in Table I, starting 
with both the position step and the momentum step, and 
the result was the same every time: If the most recently 
available version of the wave function is used whenever 
gl'^P is needed alongside the potential energy, the split- 
operator method for the GPE has the same order of accu- 
racy in the time step as the corresponding split-operator 
method for the Schrodinger equation. Whether using the 
most recent update of the wave function is also a neces- 
sary condition for the same order of time stepping error 
in the linear and in the nonlinear problem depends on 
the split-operator method on hand. For instance, if one 
starts the three-split method with a position step and 
writes 
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all choices with C2 ~ ±1 and ci — — cq give an error 
0(/i3). 

In the seven-step 0{h^) method on the order of 20000 
terms must identically cancel to validate our result, so 
clearly it is not an accident. We conjecture that our ob- 
servation holds for every minimal split-step method, for 
arbitrary high orders. Also, so far our examples have 
been in one spatial dimension. The split-operator al- 
gorithms for the Schrodinger equation work for arbitrary 
operators A and B, and go over unchanged to any number 
of spatial dimensions. Replacing the second derivative in 
position with a multidimensional Laplacian should not 
change anything in the underlying structure of the split- 



operator method or the exact result in the GPE case ei- 
ther. Our best guess is that the split-operator algorithms 
for the GPE work the same way in more than one spa- 
tial dimension. We have verified this explicitly for the 
0{h^) three-exponential splits in two and three spatial 
dimensions. 

Our observation probably originates from some alge- 
braic structure in the GPE, but at this time we can- 
not tell what structure and how. We see potential here 
for path-breaking insights into norm-preserving nonlinear 
differential equations, and maybe finite-difference equa- 
tions as well. Moreover, computationally demanding 
physics applications of the GPE arise with sets of coupled 
GPEs encountered in multi-component condensates and 
in atom-molecule coupling. We wonder if the last-update 
rule could apply quite generally to norm-preserving al- 
gorithms for such problems. One could investigate this 
in any given special case using Mathematica, but un- 
derstanding the underlying mathematics could lead to 
more penetrating answers. Integrating the GPE and re- 
lated equations in imaginary time, a frequently employed 
method to find the ground state for classical nonlinear 
fields, invites additional interesting and possibly impor- 
tant questions of the same ilk. 

In conclusion, by using the most recent existing update 
for the wave function whenever the square of the wave 
function is needed alongside the potential energy, in all of 
the cases we have studied the split-step Fourier transform 
method goes over from the Schrodinger equation to the 
Gross-Pitaevskii equation with its favorable properties 
intact. We hope that our observation is useful in future 
studies based on the Gross-Pitaevskii equation and its 
variations and generalizations. 
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